////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//---------------------------- Ebrahim Foulaadvand, 30 July 2013 -----------------                        //
//                                                                                                        //
//  The routine "JacobiNeumann" solves the 2D Laplace equation using Jacobi relaxation method.            //
//  Neumann boundary condition is used. The region is a rectangle of sides lengths Lx and Ly.             //
//                                                                                                        //
//                                                                                                        //
//                                                                                                        //                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{

const Nxsize=100,Nysize=200;

int i,j,n,Nx=100,Ny=200,T=7000;

double Lx=1,Ly=2,a1=0.,a2=0.,a3=0.1,a4=0.,delx,dely;

double Phi[Nxsize+1][Nysize+1],Phinew[Nxsize+1][Nysize+1]; // 2D arrays "Phi" and "Phinew" store the current and updated values of potential at grid points.


ofstream file1 ("Y profile x=0.5Lx n=10.plt");     // output file for the potential profile at timestep n=10.

ofstream file2 ("Y profile x=0.5Lx n=50.plt");     // output file for the potential profile at timestep n=50.

ofstream file3 ("Y profile x=0.5Lx n=300.plt");    // output file for the potential profile at timestep n=50.

ofstream file4 ("Y profile x=0.5Lx n=500.plt");    // output file for the potential profile at timestep n=50.


delx=Lx/Nx; cout<<"delx= "<<delx<<endl;

dely=Ly/Ny; cout<<"dely= "<<dely<<endl; getch();


// ----------------------- Assignment of boundary values ----------------------------------------

for(i=0;i<=Nx;i++){ Phi[i][0]=a1;  }

for(j=1;j<Ny;j++){ Phi[0][j]=a4; Phi[Nx][j]=a2;   }

Phi[0][Ny]=Phi[0][Ny-1]+a3*dely;  Phi[Nx][Ny]=Phi[Nx][Ny-1]+a3*dely;


for(i=0;i<=Nx;i++){ Phinew[i][0]=a1; }

for(j=0;j<=Ny;j++){ Phinew[0][j]=a4; Phinew[Nx][j]=a2; }

Phinew[0][Ny]=Phinew[0][Ny-1]+a3*dely;  Phinew[Nx][Ny]=Phinew[Nx][Ny-1]+a3*dely;


// ----------------------- initial potential assignement ----------------------------------------

for(j=1;j<Ny;j++){ for(i=1;i<Nx;i++){ Phi[i][j]=0.1*a3; }}


for (n=1;n<=T;n++){ if (fmod(n,100)==0) cout<<"n= "<<n<<endl; //getch();

for(j=1;j<Ny-1;j++){
for(i=1;i<Nx;i++){

Phinew[i][j]=0.25*(Phi[i+1][j] + Phi[i-1][j] + Phi[i][j+1] + Phi[i][j-1]);

}}

for(i=1;i<Nx;i++){

Phinew[i][Ny-1]=(1/3.)*(Phi[i+1][Ny-1] + Phi[i-1][Ny-1] + Phi[i][Ny-2] + a3*dely);

Phinew[i][Ny]=Phinew[i][Ny-1] + a3+dely;

}



if (n==10)    { for(j=1;j<Ny;j++){ file1<<j*dely<<" "<<Phinew[Nx/2][j]<<endl; }}

if (n==50)    { for(j=1;j<Ny;j++){ file2<<j*dely<<" "<<Phinew[Nx/2][j]<<endl; }}

if (n==300)    { for(j=1;j<Ny;j++){ file3<<j*dely<<" "<<Phinew[Nx/2][j]<<endl; }}

if (n==500)    { for(j=1;j<Ny;j++){ file4<<j*dely<<" "<<Phinew[Nx/2][j]<<endl; }}

for(j=1;j<Ny;j++){
for(i=1;i<Nx;i++){

Phi[i][j]=Phinew[i][j];

}}


}

cout<<"end." <<endl; getch();

}